Quarkonium above deconfinement as an open quantum system 
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Abstract 

Quarkonium at temperatures above deconfinement is modeled as an open quantum system, whose 
dynamics is determined not just by a potential energy and mass, but also by a drag coefficient 
which characterizes its interaction with the medium. The reduced density matrix for a heavy 
particle experiencing dissipative forces is expressed as an integral over paths in imaginary time and 
evaluated numerically. We demonstrate that dissipation could affect the Euclidean heavy-heavy 
correlators calculated in lattice simulations at temperatures just above deconfinement. 
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I. INTRODUCTION 



Because of the masses of heavy quarks relative to the QCD energy scale Aqcd, heavy- 
heavy bound states are systems which can be described to some accuracy with first quanti- 
zation techniques [IH3]. At zero temperature, lattice calculations of a heavy- heavy system 
give a potential term which can describe quarkonium spectroscopy. The potential has a non- 
trivial temperature dependence. For this reason, it has been the hope of both theorists and 
experimentalists for over twenty years that changes in the yields of quarkonium in heavy-ion 
collisions would be a signal for changes in the temperature and phase of the medium [3]. 
Since then, analysis of heavy-ion experiments have shown that some suppression of J/ip 
yields at the RHIC is anomalous [5]. However, there is little agreement among theorists 
over the cause of this suppression. It is clear that the dynamics of quarkonium above the 
deconfinement phase transition should be considered carefully. 

One step in this direction was made by Shuryak and one of us, who modeled charmo- 
nium in heavy-ion collisions as undergoing Brownian motion, with the heavy quark diffusion 
coefficient Dh taken to be quite small, as expected from phenomenology of single heavy 
quarks [9] and by gauge-gravity duality [T0HT2] . This work [7j demonstrated that the 
survival of J/ip states in plasma cannot be determined just by examining whether or not the 
temperature-dependent potentials allow bound states; instead, the dynamics of charm and 
charmonium are determined by multiple interactions with the medium, and the time scales 
of the collision are very important. 

Another important step in understanding the dynamics of charmonium is provided 
by newer lattice calculations [6]. In these calculations, the Euclidean correlators for 
local operators related to quarkonium spectroscopy are calculated, namely, G(r, T) = 
J d 3 x (J#(x, r) Jh(0, 0)), where J#(x, r) = ^(x, r)r^(x, r) is a local composite operator 
related to a given mesonic channel. Taking the results of these lattice simulations, Mocsy 
and Petreczky compared these correlators, at various temperatures and in different channels, 
with the results of two different potential models [13] . In this work, there is a low- frequency 
contribution to the spectral densities used to determine the Euclidean correlator, which is 
due to the diffusion of heavy quarks [H] . however, for u ^> T 2 /M, the diffusion of heavy 
quarks is assumed not to have any effect on the dynamics of quarkonium. 

In this paper, we will outline how the imaginary-time propagator with periodicity (3 = 1/T 
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can be determined for quarkonium systems interacting with a heat bath. In Section |TT| we 
review the reduced density matrix, apply the model of Caldeira and Leggett for quantum 
Brownian motion to imaginary time, and then describe briefly path-integral Monte Carlo 
techniques for calculating this reduced density matrix. In Section |III[ we show how the 
imaginary-time propagator for charmonium may be calculated to high accuracy with these 
techniques, for any potential or drag coefficient, and we show the results for G(r, T) for the 
vector channel at temperatures just above the deconfinement transition. 

Not many before have considered the effect of the medium on heavy quark and quarko- 
nium Euclidean correlators calculated on the lattice. Current work by Beraudo et al. [TS] 
does consider the medium effect on heavy quark correlators. Starting from the QCD La- 
grangian, they make the HTL approximation for the heavy quark's interaction with the light 
degrees of freedom at finite temperature, and then they integrate these degrees of freedom 
out. The focus of our work is on a heavy particle interacting quite strongly with the medium, 
as suggested by calculations taking advantage of gauge-gravity duality [10J. We show how 
the Lagrangian for a particle interacting strongly with a heat bath can be renormalized so 
that functional integration is well-defined, and where analytic results are possible. 

II. THE PATH INTEGRAL FORM FOR THE REDUCED DENSITY MATRIX 
OF A DISSIPATIVE SYSTEM 

For an excellent review of the functional integral approach to quantum Brownian motion, 
in both imaginary and real time, see [16]. For a briefer discussion of these ideas and how 
they apply to quarkonium, see [T7]. These results are discussed in generality, and shown to 
arise from the Schwinger-Keldysh contour integral for a heavy quark's real-time partition 
function, by Son and Teaney [18]. These systems can be studied with approaches besides 
the path-integral formulation, and the study of quantum Brownian motion, in terms of a 
partial differential equation for the density matrix, has been studied by Hu, Paz, and Zhang 

[m [20]. 
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A. The reduced density matrix 



Without loss of generality, consider a system consisting of a heavy particle of mass M 
which we will call the system, minimally coupled to a harmonic oscillator of mass m. 

L = Ls + Lr: 



L s = -Mx 2 -V(x), 

Lj = -mr 2 mu 2 r 2 — Cxr. 

2 2 



We analytically continue this Lagrangian to r = it. 



(1) 



Sg[x] = 

Sf[x,r] = I 
Jo 



-Mx 2 + V(x) 



-mr 2 H — muj 2 r 2 + Cxr 
2 2 



dr 



(2) 



We can simplify this Lagrangian by making a change of variables. We subtract the particular 
solution to the classical equations of motion as determined by Eq. [2j 

C 



r(r) = r'(r) H / dr'x(r') sinh [u(t — r')l 

mu J 

r'(r) + A[x,T\ 



(3) 



In terms of the shifted coordinate r' the Euclidean action becomes that of a simple harmonic 
oscillator: 



Sf[x,r] 



111 

-mr' 2 + -mu 2 r' 2 + -Cx(r)A[x, r] 



dr 



+mA[x,P] lr'(P) + -A[x,P] 



(4) 



As always, the propagator of a system for imaginary time f3 = —i/T gives matrix elements 
of the thermal density operator. In our example, the density matrix has 4 indices, 2 for each 
degree of freedom of the system. In the example given above the density matrix is given as 



p(x i ,r i ;x f ,r f ;f3)= I Vxl Vrexp (S§[x] - Sf [x,r]) . 

'x(0)=Xi Jr(0)=ri 



(5) 



If we were never interested in measurements of the degree of freedom r, we could take the 
trace over the indices corresponding to this degree of freedom, and work with a density 



operator with only two indices. With this in mind, we define the reduced density matrix: 

p red (xi,x f ,/3) = / drp(xi,r;x f ,r; (3). (6) 



In the next section it will be the only operator practical for calculating thermal averages. 

For the system defined by Eq. [TJ we can now write a path-integral description of the 
reduced density matrix, 

p red (xi,x f ,f3) = J dr J Vx Vr' exp (S§[x] - Sf [x,r')) 

= J Vxexp(-SZ{ x }) J dr J Vr'exp^-S'fix,^ , (7) 

where we integrate over paths with endpoints x(0) = x i: x(j3) = Xf, r'(0) = r, and r'(/3) = 
r — A[x,/3]. The final step in [7] allows the integral over the paths r'(r) to be performed 
independently of the integral over x(r). Using Eq. |1J a simple Gaussian integral yields 

rx(/3)=Xf 

p red (xi,Xf,f3) = Vx 

J x(0)=Xi 

x e xp(-Sf + °\, uli ^ f dr [ ds x(r)x(s) cosh [c* (r - s - /3/2)] ) (8) 

where the summation is introduced because we have generalized this to a system where a 
heavy particle interacts with a bath of simple harmonic oscillators. This is now the path 
integral form for the reduced density matrix. This is analogous to the idea of influence 
functionals, worked out for real-time propagators many years ago [2Tj . 

B. Making a dissipative system 

Any finite quantum-mechanical system is reversible and therefore inappropriate for de- 
scribing Brownian motion. One might find it intuitive that if the bath of harmonic oscillators 
were taken to an infinite limit, it would be "large enough" so that energy from the heavy 
particle could dissipate into the system and never return. This intuition was proven to be 
true by the authors of |22j, who considered the real-time evolution of the density matrix 
for our system and showed that when the bath of harmonic oscillators is determined by the 
continuous density of states 

C 2 (u)p D (to) = { * ~ (9) 

if to > n 



the force autocorrelator for the heavy particle is proportional to 5{t—t') at high temperatures. 
In this "white noise" limit, the density matrix describes an ensemble of particles interacting 
according to the Langevin equation, which has been used to describe Brownian motion for 
a long time. It is a stochastic differential equation, is irreversible, and evolves any ensemble 
towards the thermal phase space distribution. 

We now take this density of states and substitute it into Eq. [8] 



p red (xi,x f ,/3) 



x exp 



S§[x] 



Vx 





x{0)=Xi 
Q 

dcu I dr I ds x(r)x(s] 
< i Jo Jo 



lo cosh [oj(t — s — (3/2)] 
sinh(f ) 



(10) 



The divergences of this action can be isolated by integrating by parts twice 
io Jo Jo 



sinh(f ; 



1 



dr (x(r)) 2 --(x i -x f ) 2 \n(Mn/ V ) 



o 



1 



(Xi 



2 

1e + In 



cosh(Q/3/2) - 
sinh(fi^/2) 



irM 



+ {Xi — Xf) / dr x(t) lnsin 
Jo 

+ I dr ds x(t)x(s) lnsin 
Jo Jo 



71T 

7r(r — s) 



(11) 



where lnsin(x) = In [sin (x)} and 7# is the Euler-Mascheroni constant. The first two terms 
on the right-hand side correspond to a renormalization of the potential for the heavy par- 
ticle, always necessary when considering the interaction of a particle with infinitely many 
additional degrees of freedom. They are temperature-independent in the limit of large Q, 
and may be introduced as temperature-independent modifications of the Lagrangian for the 
particle. The final three terms are finite and well-behaved, and have been evaluated in the 
limit Q — > oo. The final form for the reduced density matrix becomes 

-iE\ 



p red (xi,Xf,/3) 



x(0)=Xi 



Vx exp ^ — S s [x 
V 



2tt 
V 



(Xi - x f f 

Xi X j 



Ie + In 



dr x(t) lnsin 



o 



+ — / dr ds x(t)x(s) lnsin 

Jo Jo 



TIT 

J 
7t(t 



(12) 



In summary, we have determined an expression for the reduced density matrix of a system 
(consisting of a massive particle in a potential) interacting with a bath of oscillators. The 
coupling to the bath has been chosen in order to reproduce the results of classical Brownian 
motion in the high temperature limit. We have taken care to make the term in the expo- 
nential finite, by isolating the divergences through integration by parts. This is important 
for path integral Monte Carlo simulation to be possible for this functional integral. 



C. Example: The otherwise free particle 



The reduced density matrix for a particle interacting with such a bath can be determined 
analytically for the otherwise free particle (Vr(x) = 0). First, write an arbitrary path as an 
expansion around the classical solution 

x(t) = x d (r) +£(t); 
x c i{r) = Xi + (x f - Xi)r/(3, 

e(r) ss £c„sin(^). (13) 



n=l 



V P 



One can find an analytic solution by substituting this into Eq. 12 . We skip the details: 



after evaluating numerous integrals using contour integration, and then changing variables 
for the integration over the even Fourier coefficients, we find the reduced density matrix 



Predy^ii ^ f i P) 



M 7] 
27T/3 + 2TX 2 



x exp < — 



M T] 

2/3 + 27 



log(2) + 7 £ + ( 1 4 

V (3 



log 



2vrM 



2vrM 
*( 1 + 



2ttM 



(xt-xtf ,(14) 



where *ff(x) is the digamma function and the overall normalization is determined by analytic 
continuation of /3 and requiring the propagator to conserve probability for purely imaginary 
p. One interesting physical result is easily obtained from the Fourier transform of the 
reduced density matrix 



(P 2 ) 



M, 



eff 



(3 



M f 



ff 



M + 



71 



log 



VP 
2nM 



V 1 



V_P \ 
2nM 



(15) 



D. The path-integral Monte Carlo algorithm for the reduced density matrix 



Obtaining the analytic result for our reduced density matrix was reasonable for the oth- 
erwise free particle. For the simple harmonic oscillator, the analytic result exists but has a 



rather complicated expression. For the potentials which describe quarkonium spectroscopy 
with some precision, analytic work becomes entirely impractical. We would like to use nu- 
merical simulation to obtain reliable estimates of reduced density matrices and Euclidean 
correlators. 

The most natural numerical approach for the formalism we adopted is path-integral Monte 
Carlo (PIMC). For an excellent review of the technique, see the review of D. M. Ceperley 
[23] . Paths are either sampled according to a Metropolis algorithm determined by the action 
of interest, or sampled from some convenient distribution which samples the entire space 
of paths with some weight. For our work with a single degree of freedom, we found that 
sampling a convenient distribution (in our case, the distribution of paths determined by 
exp(— Sf ree [x})) to be sufficient, which is easily sampled for any discretization of the path 
with a bisection method. When sampling the space of paths, the next step is to determine 
an estimate for the action of each path. For our case, the primitive action with the simplest 
integration of the new dissipative terms is sufficient. Once this is determined, any correlator 
can be calculated by sampling paths and making weighted averages. 

III. QUARKONIUM AS AN OPEN QUANTUM SYSTEM 

We should now apply this functional integral to the problem of quarkonium at tempera- 
tures above deconfinement. To this end, let us rework the Euclidean correlators calculated 
on the lattice to a form with which we can calculate. 

For a given channel, the two-point Euclidean correlator for a composite mesonic operator 
is given by 



where f2 = 1, 7 , 7^, 7°7^ determine the mesonic channel to be scalar, pseudo-scalar, vector, 
or pseudo-vector, respectively. 

For now, consider only the vector channel (the following arguments must be modified for 
the scalar and pseudo- vector channels). In order to switch to a first-quantization approach 
for the energy range where a potential model would be appropriate, think of this correlator as 
being the sum of the expectation values of an operator over all of the states in the Fock space 





Jn(x,r) = ^(x, r)fh/>(x,r) 



(16) 
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of iV-particle mesonic systems, with each expectation value entering the sum weighted by 
the state's Boltzmann factor. One term in this sum, of course, is the "vacuum" expectation 
value 

(0|J"(x,r)^(O,0)|0> = (//,x;T|//,O;0) light , (17) 

where the subscript represents that the trace is still taken over all of the light degrees of 
freedom in QCD. 



The right-hand side of Equation [17] represents an important step: by decomposing the 
mesonic operators into products of creation and annihilation operators, the vacuum expecta- 
tion value can be related to the overlap of two one-particle states. Also, note here that since 
expectation values of the higher states enter into the thermal average multiplied by factors 
which are roughly powers of exp(— 2Mq/T), the vacuum expectation value dominates the 
thermal average in the limit M ^> T. Therefore, we have identified the leading contribution 
to the correlator in the vector channel. 

The trace over the light degrees of freedom is of course the central problem of QCD, and 
no analytic result exists. However, in the infinite mass limit, this has been computed on the 
lattice as the expectation value of two Polyakov loops [21]. The dissipative effects on this 
propagator have been studied, as we noted previously, by gauge-gravity duality. Our ansatz 
here is that these results describe well the dynamics of sufficiently heavy quarks. 

The path integral we wish to evaluate can be written in terms of relative and absolute 
coordinates as 

7T 



(x,x;r|O,O;0) = / DX exp 
x / Ux re i exp 



7T 



dr 



ds X(r)X(s) lnsin 







r - s 



+ 



2tt 







dr 
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f 


dr [ 


■Jo 


Jo 



ds ± re i(T)± rel (s) lnsin 



7T 



where the functional integration are over all paths satisfying X(0) = 0, X(r) = x, X(/3) = 0, 
x rd (0) = 0, x reZ (r) = 0, and x rei (/3) = 0. 



The right-hand side of Equation 18 can now be calculated, for any potential, with a PIMC 
simulation. We use the Cornell potential for the interaction between the heavy quarks. 



To obtain G(r), the first term on the right-hand side of Eq. 18 is integrated with respect 
to X, yielding the diagonal value for the free-particle reduced density matrix which is inde- 
pendent of r. In Figure [T] we show the mesonic correlator as a function of r at T = 1.2T C . 
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FIG. 1: G(t)/G v=q evaluated for rj = 0.1 GeV~ 2 at T = 1.2T C . 

We have chosen a representative value of 7/ = 0.1 GeV~ 2 and show results for fixed bare and 
effective charm mass. 

Clearly, dissipative effects due to the medium change the mesonic correlators measured 
in lattice simulations. Also, these effects are not trivial: holding the effective charm mass 
fixed but varying r\ still leads to changes in the mesonic correlators. The next step is to 
perform spectral analysis of these correlators to examine how dissipation affects the spectral 
function for this correlator. From here, the dissipative effects can be included in heavy-ion 
phenomenology. 

IV. CONCLUSIONS 

We have taken the approach of Caldeira and Leggett and determined an imaginary-time 
path integral, representing the reduced density matrix for a dissipative system, which con- 
tains integrals that are convergent and therefore suitable for numerical integration. We have 
shown how interactions with the QCD medium above deconfinement can affect Euclidean 
heavy-heavy correlators. 

The purpose of this paper is only to identify this effect. In order to achieve a result which 
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is appropriate for comparisons with the lattice data, the following work must be done: 1.) 
the mass of the heavy quark must be changed whenever the drag coefficient is changed, since 
we are concerned with results for G(t) with different values of rj but the same "effective" 
mass, 2.) the maximum entropy method must be applied to G(t), so that the spectral 
function corresponding to the correlation function will be extracted with high precision, 3.) 
this spectral function must have the zero-mode contribution added and be cut so that the 
perturbative result for the spectral function is used at large u instead, and finally 4.) this 
spectral function should be used to calculate G(t). However, let us note here that this work 
can be put to greater use in examining quarkonium in heavy-ion collisions. 

Let us argue one last time for treating quarkonium above deconfinement as an open 
quantum system: typically in particle physics, the rate for the scattering of an iV-particle 
state into another is determined from the square of the matrix element whose indices are 
the initial and final states, where these states are in the momentum basis. This makes per- 
fect sense in high-energy experiments, where the incoming and outgoing states are basically 
each 2-particle momentum eigenstates and the matrix element can be expanded in terms 
of a small coupling. Such an approach, however, seems entirely inappropriate for quarko- 
nium formed in a heavy-ion collision, whose constituent quarks are localized in position 
relative to the surrounding medium and are strongly coupled both to the medium and each 
other. Non-perturbative techniques and experiment both suggest that quarkonium rapidly 
thermalizes and interacts strongly with the medium. The most reasonable approach for 
explaining the observables is to describe quarkonium with a reduced density matrix, whose 
evolution is determined by a potential and a drag coefficient which are both determined 
non-p ert urb at ively. 
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